#!/home/timm/opt/gawk/bin/gawk -f

BEGIN {
    TooMuch  = 20
    Enough   = 3
    Sampleth = 10
    Prune    = 1
    FS=","; OFS=",\t"
    main()
}
function main(     w,gp,np,hi,lo,rows,all,d) {
    while(getline > 0) {
	if (NR==1) {
	    srand(Seed ? Seed : 1); 
            header(w,gp,np,hi,lo,rows,all,d)
	} else 
	  rows=row(w,gp,np,hi,lo,rows,all,d) 
    }
    think(w,gp,np,hi,lo,rows,all,d)
}
function header(w,gp,np,hi,lo,rows,all,d,   i) {
    for(i=1;i<=NF;i++)
	if ($i !~ /\~/) {
	    w[i]  = 1
	    gp[i] = $i ~ /!/
	    np[i] = $i ~ /\$/
	    if (np[i])
		hi[i] = -1* (lo[i]=10^32)}
}
function row(w,gp,np,hi,lo,rows,all,d,   i) {
    rows++
    all[rows]=rows
    for(i in np) {
	d[rows][i]=$i
	if ($i != "?")
	    if (np[i]) {
		if ($i > hi[i]) hi[i] = $i
		if ($i < lo[i]) lo[i] = $i }}
    return rows
}
function think(w,gp,np,hi,lo,rows,all,d,\
	       min,more,rounds,err,id,r,x) {
    min  = 10^32
    more = Enough
    id   = 1
    while(more) {
        if (++rounds > TooMuch) break
	say(".")
	err = project(x,id,min,w,gp,np,hi,lo,rows,all,d)
	if (err < min) {
	    more--
            id++
	    print (min=err) " rnn" 
	} else
	    delete x[id]
    }
    print ""
    for(id in x) rnn(id,x[id],r)
    for(id in r) xx[r[id]]++
    for(id in xx)  {
	per = int(100*xx[id]/length(all))
	print id,xx[id],per "%",stars(per) | "sort -t, -n "
    }

}
function xsort(r1,x1,r2,x2) {
    return x1["x"] - x2["x"]
}

function rnn(id,a,r,    i,j,k,b,n,right,left) {
    n=asort(a,b,"xsort") 
    r[2]++
    r[n-1]++
    for(j=2;j<n;j++) {
       i = j-1
       k = j+1
       r[k] += 0; r[i] += 0; r[j] += 0
       left = b[j]["x"] - b[i]["x"]
       right= b[k]["x"] - b[j]["x"]
       if      (left > right)  { r[k]++ }
       else if (left < right)  { r[i]++ }
       else                    { r[k]++; r[i]++ }
    }
    # print ""
    # for(j  =1;j<=n;j++) {
    # 	com = "sort -t, -n -k 4 -n -k 3"
    # 	at = int(b[j]["x"]*100)
    # 	print id,at, b[j]["id"],r[j],stars(r[j])  | com
    # }
    # close(com)

}

function project(xs,id,min,w,gp,np,hi,lo,rows,all,d,	\
                 one,two,i,x,y,a,b,c,row,err,m,n) {
    m    = length(all)
    one  = any(all)
    do { two = any(all) } while (two == one)
    c = dist(one,two,w,gp,np,hi,lo,rows,all,d)
    for(i in all) {
        row = all[i]
	if (row == one ) {
	    a=x=0; b = c
        } else {
	    if (row == two) { 
		a=x=c;  b = 0
	    } else {
		a = dist(row,one,w,gp,np,hi,lo,rows,all,d)
		b = dist(row,two,w,gp,np,hi,lo,rows,all,d)
		x = (a^2+c^2 - b^2)/(2*c + 0.00001);  
	    }}	   
	y = sqrt(abs(a^2 - x^2))
        err += y
        xs[id][i]["x" ] = x
	xs[id][i]["id"] = i
        if (Prune) {
	    n++
	    if (n > Sampleth && n > m/Sampleth)
		if (err/n > min) 
		    return 10^32 } 
    }
    return err/m
}
function dist(row1,row2,w,gp,np,hi,lo,rows,all,d,\
              col,sum,ws,w1,x,y) {
    if (row1==row2) 
	return 0
    for(col in np) {
	if (gp[col]) continue
	w1 = w[col]; 
	if (w1==0) continue
	x  = d[row1][col]
	y  = d[row2][col]
	if (x=="?" && y=="?") continue 
	sum += w1 * dist1(x,y,col,np,hi,lo)
	ws  += w1
    }
    return sqrt(sum)/sqrt(ws + 0.000001)
}
function dist1(x,y,col,np,hi,lo) {
    if (x == y) 
	return 0 
    if (!np[col])  
	return value(x) != value(y)
    if (y == "?") {
	x  = norm(x,col,hi,lo)
	y  = x > 0.5 ? 0 : 1
    } else if (x == "?") {
	y = norm(y,col,hi,lo)
	x = y > 0.5 ? 0 : 1
    } else {
	x = norm(x,col,hi,lo)
	y = norm(y,col,hi,lo)
    }
    return  (x - y)^2 
}
function norm(x,col,hi,lo) {
  return (x - lo[col])/(hi[col] - lo[col] +0.00001)
}
function value(x) {
    return x=="?" ? "SomeCrazySyMBOL" : x
}
function abs(x) { return x < 0 ? -1*x : x }
function any(l, n,max) {
    max = length(l)
    n   = int(rand()*max+0.5)
    if   (n < 1)      return l[1]
    else if (n > max) return l[max]
    else return l[n]
}
function say(x) {
    printf(x)>"/dev/stderr"
    fflush("/dev/stderr")
}
function o(a,s,  i){
    for(i in a) print s"[" i "]=" a[i]
}
function l2sd(a,    i,n,x,delta,mean,m2) {
    for(i in a) {
	n    += 1
	x     = a[i]
	delta = x - mean
	mean += delta/n
	m2   += delta*(x -mean)
    }
    return sqrt(m2/(n-1))
}
function stars(n,s,   out) {
    s = s ? s : "*"
    while (n-- > 0) out = out s
    return out
}